﻿using System;
using System.Collections.Generic;
using System.Text;
using System.Windows.Forms;

namespace GeoFly
{
/// <summary>
/// 栅格产流计算时用的类
/// </summary>
public class GridWaterBalance
{

    public GridWaterBalance()
    {
    }
	
	/// <summary>
	/// 设置栅格参数
	/// </summary>
	/// <param name="currow"></param>
	/// <param name="curcol"></param>
	/// <param name="rint">降雨强度</param>
	/// <param name="dFp">土壤下渗率</param>
	/// <param name="curDate"></param>
	/// <param name="hr"></param>
    public void SetGridPara(int currow,int curcol,double rint,double dFp, DateTime curDate, double hr)
    {
        this.row = currow;
        this.col = curcol;
        this.m_dRIntensity = rint;
        this.m_nYear = curDate.Year;
        this.m_nDn = DateAndTime.GetDnInYear(curDate.Year,curDate.Month,curDate.Day);
        this.m_dFp = dFp;
        this.m_dHr = hr;

        this.m_dPcp =HydroSimulate.g_ClimatePara.pcpdata.m_gridLayer.Values[row, col];
        this.m_dHeight =HydroSimulate. g_GridLayerPara.g_DemLayer.Values[row, col];
        int SoilID = (int)HydroSimulate.g_GridLayerPara.g_SoilLayer.Values[row, col];
        this.m_dFc = HydroSimulate.g_GridLayerPara.SoilTypes[SoilID].SP_Stable_Fc;

        if (HydroSimulate.g_ModelRunPara.RunoffSimuType ==RunOffSimuType.STORM_RUNOFF_SIMULATION)
            return;

        if (HydroSimulate.g_ClimatePara.PetInfo.DataFileName == "")
            this.m_dTav = HydroSimulate.g_ClimatePara.tmpmeandata.m_gridLayer.Values[row, col];

        this.m_nMon =DateAndTime.GetMonthByDn(m_nYear, m_nDn);
        this.m_dGridLAI = this.DailyLAI(HydroSimulate.g_ModelRunPara.DLAICalcMethod,HydroSimulate.g_GridLayerPara);
        this.m_dGridCovDeg = this.DailyCoverDeg(HydroSimulate.g_GridLayerPara);
    }
	/// <summary>
	/// 计算冠层截留
	/// </summary>
	/// <param name="g_ModelRunPara"></param>
	/// <param name="g_ClimatePara"></param>
    public double CalcCI()
    {
        m_CanopyStore.SetGridValue(m_dGridLAI, m_dPcp, m_dGridCovDeg, 0.046);
        m_dCrownInterc = m_CanopyStore.CanopyStore();
        return m_dCrownInterc;
    }
	/// <summary>
	/// 计算潜在蒸散发
	/// </summary>
	/// <param name="g_ModelRunPara"></param>
	/// <param name="g_ClimatePara"></param>
	/// <param name="dalbedo">默认值为0.23</param>
    public double CalcPET(ModelRunPara g_ModelRunPara, ClimatePara g_ClimatePara, double dalbedo)//dalbedo的默认值为0.23
    {
        m_dPET = 0;
        if (g_ModelRunPara.PETMethod == PET_METHOD.PET_REAL)
        {
            m_dPET = g_ClimatePara.petdata.m_gridLayer.Values[row, col];
        }
        //如果最高气温和最低气温的文件不存在，则退出执行该函数
        if (g_ClimatePara.TempmxInfo.DataFileName == "" || g_ClimatePara.TempmnInfo.DataFileName == "")
            return -9999;
        double dTmx = g_ClimatePara.tmpmxdata.m_gridLayer.Values[row, col];//最高气温
        double dTmn = g_ClimatePara.tmpmndata.m_gridLayer.Values[row, col];//最低气温


        switch (g_ModelRunPara.PETMethod)
        {
            case PET_METHOD.PET_PRISTLEY_TAYLOR:
                {
                    PETInPristley prist = new PETInPristley(m_dTav, m_dHeight, m_nYear, m_nDn, 33);
                    double dRLong = prist.NetLongWaveRadiationRHmd(g_ClimatePara.slrdata.m_gridLayer.Values[row, col],
                                                            g_ClimatePara.hmddata.m_gridLayer.Values[row, col]);
                    double dRShort = prist.NetShortWaveRadiation(dalbedo, g_ClimatePara.slrdata.m_gridLayer.Values[row, col]);
                    m_dPET = prist.PETByPristley(1.26, 0);
                }
                break;
            case PET_METHOD.PET_HARGREAVES:
                {
                    PETInHargreaves har = new PETInHargreaves(m_dTav, m_dHeight, dTmx, dTmn, m_nYear, m_nDn);
                    double dRLong = har.NetLongWaveRadiationRHmd(g_ClimatePara.slrdata.m_gridLayer.Values[row, col],
                                                            g_ClimatePara.hmddata.m_gridLayer.Values[row, col]);
                    double dRShort = har.NetShortWaveRadiation(dalbedo, g_ClimatePara.slrdata.m_gridLayer.Values[row, col]);
                    m_dPET = har.PETByHarg();
                }
                break;
            case PET_METHOD.PET_FAO_PENMAN_MONTEITH:
                {
                    PETInFAOPM faopm = new PETInFAOPM(m_dTav, m_dHeight, dTmx, dTmn, m_nYear, m_nDn, 33);
                    double dRLong = faopm.NetLongWaveRadiationRHmd(g_ClimatePara.slrdata.m_gridLayer.Values[row, col],
                                                            g_ClimatePara.hmddata.m_gridLayer.Values[row, col]);
                    double dRShort = faopm.NetShortWaveRadiation(dalbedo, g_ClimatePara.slrdata.m_gridLayer.Values[row, col]);
                    double G = 0;
                    m_dPET = faopm.PETByRAVP(g_ClimatePara.WindData.m_gridLayer.Values[row, col],
                                            g_ClimatePara.hmddata.m_gridLayer.Values[row, col], G);
                }
                break;
            case PET_METHOD.PET_DEBRUIN:
                {
                    PETInDeBruin debruin = new PETInDeBruin(m_dTav, m_dHeight, m_nYear, m_nDn, 33);
                    double dRLong = debruin.NetLongWaveRadiationRHmd(g_ClimatePara.slrdata.m_gridLayer.Values[row, col],
                                                                g_ClimatePara.hmddata.m_gridLayer.Values[row, col]);
                    double dRShort = debruin.NetShortWaveRadiation(dalbedo, g_ClimatePara.slrdata.m_gridLayer.Values[row, col]);
                    m_dPET = debruin.PETByDeBruin();
                }
                break;
            case PET_METHOD.PET_REAL:
                {
                    m_dPET = g_ClimatePara.petdata.m_gridLayer.Values[row, col];
                }
                break;
            case PET_METHOD.PET_PENMAN_MONTEITH:
                {
                    PETInPM pm = new PETInPM(m_dTav, m_dHeight, m_nYear, m_nDn, 33);
                    double dRLong = pm.NetLongWaveRadiationRHmd(g_ClimatePara.slrdata.m_gridLayer.Values[row, col],
                                                        g_ClimatePara.hmddata.m_gridLayer.Values[row, col]);
                    double dRShort = pm.NetShortWaveRadiation(dalbedo, g_ClimatePara.slrdata.m_gridLayer.Values[row, col]);
                    m_dPET = pm.PETInPMByRHmd(g_ClimatePara.WindData.m_gridLayer.Values[row, col],
                        g_ClimatePara.hmddata.m_gridLayer.Values[row, col], 0);
                }
                break;
        }
        return m_dPET;
    }
														
	/// <summary>
	/// 计算实际蒸散发
    /// 由于Kojima法会出现比较大的空间不连续性，所以
    /// 目前只有互补相关理论法可以用。	
	/// </summary>
	/// <param name="g_ModelRunPara"></param>
	/// <param name="g_ClimatePara"></param>
	/// <param name="dalbedo"></param>

    public void CalcAET(ModelRunPara g_ModelRunPara, ClimatePara g_ClimatePara, double dalbedo)//dalbedo的默认值为0.23
    {
        m_dAET = 0;
        if (m_dPET == 0) return;

        switch (g_ModelRunPara.AETMethod)
        {
            case AET_METHOD.AET_BY_CROP_COEFFICIENTS:
            case AET_METHOD.AET_BY_COMPRELATIONSHIP:
                {
                    m_dAET = CompleAET(dalbedo);
                    break;
                }
            case AET_METHOD.AET_BY_COMPRELA_AND_KOJIMA:
                {
                    double dtempmean = g_ClimatePara.tmpmeandata.m_gridLayer.Values[row, col];

                    if (dtempmean < 0)
                    {
                        WaterVapor waterv = new WaterVapor(dtempmean);
                        double dSatVap;
                        double dActVap;
                        dSatVap = waterv.SatuVapPressure();
                        dActVap = waterv.ActVapPressure(g_ClimatePara.hmddata.m_gridLayer.Values[row, col]);
                        m_dAET = 0.001 * g_ClimatePara.WindData.m_gridLayer.Values[row, col]
                                    * (dSatVap - dActVap) * 10 * 240;
                    }
                    else
                        m_dAET = CompleAET(dalbedo);
                }
                break;
            default:
                break;
        }
    }
													
	/// <summary>
	/// 计算干旱指数
	/// </summary>
	public void	CalcAI()
    {
        m_dAI = 0;
    }
	/// <summary>
	/// 计算净雨(计算冠层截留量与降雨量之间的差)，若冠层截留量小于降雨量，为负值
    /// 否则，为正值
	/// </summary>
	public void	CalcNetRain()
    {
        this.m_dNetRain = m_dPcp;
	    if(m_dPcp < m_dCrownInterc)
	    {
		    m_dCIDeficit = m_dCrownInterc - m_dPcp;
		    m_dNetRain = 0;
	    }
	    else
	    {
            m_dCIDeficit = 0;
		    m_dNetRain = m_dPcp - m_dCrownInterc;
	    }
    }
    /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
    +														+
    +				功能：计算栅格总径流					+
    +		(原理见博士论文：栅格水量平衡关系计算产流)		+
    +														+
    +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    +														+
    +		总径流 ＝ 地表径流 ＋ 壤中流 ＋ 地下径流		+
    +			1、计算地表径流								+
    +			2、计算壤中流								+
    +			3、计算地下径流								+
    +			4、计算土壤含水量变化						+
    +														+
    +		返回值：栅格产流类型；							+
    +														+
    +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    +														+
    +		调用本函数前的主要输入参数：					+
    +			m_dNetRain	  -- 栅格净雨；					+
    +			m_dAET		  -- 实际蒸散量；				+
    +			m_dRIntensity -- 雨强；						+
    +			m_dFc		  -- 土壤稳定下渗率；			+
    +			m_dFp		  -- 时段土壤下渗率；			+
    +														+
    +++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
	/// <summary>
	/// 计算产流类型，并同时生成地表、壤中、地下径流和基流
	/// </summary>
	/// <param name="g_GridLayerPara"></param>
	/// <param name="g_ModelInputPara"></param>
    /// /// <param name="SP_Sw">栅格的当前土壤含水量</param>
	/// <returns></returns>
    public int CalcRunoffElement(GridLayerPara g_GridLayerPara, ModelInputPara g_ModelInputPara,double SP_Sw)
    {
        int waterYieldType = 0;  
        //获取土壤水赤字
        int SoilID=(int)g_GridLayerPara.g_SoilLayer.Values[row, col];
        //包气带土壤缺水情况
        double dthet = g_GridLayerPara.SoilTypes[SoilID].SoilWaterDeficitContent(SP_Sw);
        double dPE = this.m_dNetRain - m_dAET;

        double dIFc = m_dRIntensity - m_dFc;//雨强减去稳定下渗量
        double dIFp = m_dRIntensity - m_dFp;//雨强减去土壤实际下渗率
        
        double dBaseQ=0;
        double dLatQ=0;
        m_dTotalQ = m_dSurfQ = m_dLateralQ = m_dBaseQ = dBaseQ = dLatQ = 0;

        if (dthet < 0) //前期土壤未蓄满
        {
            if (this.m_dNetRain < this.m_dAET) //净雨量小于蒸发量
            {
                if (dIFc < 0)//不超渗
                {
                    waterYieldType = 1;//方案1：不产流
                }
                else  //超渗								
                {
                    m_dSurfQ = dIFc * m_dHr;
                    waterYieldType = 2;//方案2：超渗产流
                }
            }
            else  //净雨量大于蒸发量
            {
                if (this.m_dRIntensity < this.m_dFc) //不超渗				
                {
                    m_dBaseQ = this.m_dFc * m_dHr;  //计算基流量
                    m_dSurfQ = dPE - m_dBaseQ; //超渗量减去基流就是表面流
                    if (m_dSurfQ < 0)
                    {
                        m_dSurfQ = dPE * g_ModelInputPara.SurfQOutFactor;
                        m_dBaseQ = dPE * (1 - g_ModelInputPara.SurfQOutFactor);
                    }
                    waterYieldType = 3;//方案3：蓄满产流
                }
                else  //超渗
                {
                    m_dSurfQ = dIFc * m_dHr;
                    m_dBaseQ = this.m_dFc * m_dHr;
                    if (m_dSurfQ + m_dBaseQ > dPE)
                    {
                        m_dBaseQ = dPE - m_dSurfQ;
                        if (m_dBaseQ < 0)
                        {
                            m_dSurfQ = dPE * g_ModelInputPara.SurfQOutFactor;
                            m_dBaseQ = dPE * (1 - g_ModelInputPara.SurfQOutFactor);
                        }
                    }
                    else
                    {
                        m_dLateralQ = dPE - m_dSurfQ - m_dBaseQ;
                    }
                    waterYieldType = 4;//方案4：超渗产流＋蓄满产流
                }
            }
        }
        else  //前期土壤已蓄满
        {
            if (dPE < dthet)   
            {
                if (m_dRIntensity < m_dFp)						
                {
                    waterYieldType = 5;    //方案5：不产流
                }
                else								
                {
                    m_dSurfQ = dIFp * m_dHr;
                    waterYieldType = 6;    //方案6：超渗产流
                }
            }
            else
            {
                if (m_dRIntensity < m_dFp)			//方案7：蓄满产流
                {
                    m_dBaseQ = m_dFp * m_dHr; //基流量就是下渗量
                    m_dSurfQ = dPE - dthet - m_dBaseQ;
                    if (m_dSurfQ < 0)
                    {
                        m_dSurfQ = (dPE - dthet) * g_ModelInputPara.SurfQOutFactor;
                        m_dBaseQ = (dPE - dthet) * (1 - g_ModelInputPara.SurfQOutFactor);
                    }
                    waterYieldType = 7;
                }
                else								
                {
                    m_dTotalQ = dPE - dthet;
                    m_dSurfQ = dIFp * m_dHr;
                    m_dBaseQ = m_dFp * m_dHr;
                    if (m_dSurfQ > m_dTotalQ)
                    {
                        m_dBaseQ = 0;
                        m_dSurfQ = m_dTotalQ;
                        m_dLateralQ = 0;
                    }
                    else if (m_dSurfQ + m_dBaseQ > m_dTotalQ)
                    {
                        m_dBaseQ = m_dTotalQ - m_dSurfQ;
                        if (m_dBaseQ < 0)
                        {
                            m_dSurfQ = (dPE - dthet) * g_ModelInputPara.SurfQOutFactor;
                            m_dBaseQ = (dPE - dthet) * (1 - g_ModelInputPara.SurfQOutFactor);
                        }
                        m_dLateralQ = 0;
                    }
                    else
                    {
                        m_dLateralQ = m_dTotalQ - m_dSurfQ - m_dBaseQ;
                    }
                    waterYieldType = 8;  //方案8：超渗产流＋蓄满产流
                }
            }
        }
        //蒸发剩余量-基流量-表面流量-壤中流量 
        this.Cell_SP_Sw = SP_Sw + dPE - m_dBaseQ - m_dSurfQ - m_dLateralQ;
        
        
        //如果土壤含水量超过调萎系数
        if(this.Cell_SP_Sw>g_GridLayerPara.SoilTypes[SoilID].SP_Wp)        
        {
            if (waterYieldType == 3 || waterYieldType == 4 || waterYieldType == 7 || waterYieldType == 8)
            {
                double dPerco = SWPercolation(this.Cell_SP_Sw,g_GridLayerPara.SoilTypes[SoilID].SP_Fc,m_dHr,g_GridLayerPara.SoilTypes[SoilID].TPercolation);

                if (dPerco > 0)
                {
                    this.Cell_SP_Sw -= dPerco;
                    m_dBaseQ += dPerco;

                    dLatQ = dPerco * g_ModelInputPara.LatQOutFactor;
                    this.m_dLateralQ += dLatQ;
                    this.Cell_SP_Sw -= dLatQ;
                }
            }
            if(this.Cell_SP_Sw>g_GridLayerPara.SoilTypes[SoilID].SP_Wp)
            {
                double dRecharge = SWRecharge(this.Cell_SP_Sw, m_dHr, m_dFp, g_GridLayerPara.SoilTypes[SoilID].SP_Init_F0);
                if (dRecharge > this.Cell_SP_Sw - g_GridLayerPara.SoilTypes[SoilID].SP_Wp)
                    dRecharge = this.Cell_SP_Sw - g_GridLayerPara.SoilTypes[SoilID].SP_Wp;

                dBaseQ = dRecharge * (1 - Math.Exp(-1.0 * (dRecharge / this.Cell_SP_Sw)));
                if (dBaseQ < 0)
                    dBaseQ = 0;

                m_dBaseQ += dBaseQ;
                this.Cell_SP_Sw -= dBaseQ;
            }

        }
        else
            this.Cell_SP_Sw = g_GridLayerPara.SoilTypes[SoilID].SP_Wp * 1.05;
        if (Cell_SP_Sw < 0)
            throw new Exception("土壤含水量不能为负值");
        double spsw = this.Cell_SP_Sw;
        if (m_dBaseQ < 0)
            MessageBox.Show("计算的地下径流分量为负。", "GridWaterBalance");
        if (m_dSurfQ < 0)
            MessageBox.Show("计算的地表径流分量为负。", "GridWaterBalance");
        if (m_dLateralQ < 0)
            MessageBox.Show("计算的壤中径流分量为负。", "GridWaterBalance");

        return waterYieldType;
    }

    double CompleAET(double dalbedo)//dalbedo的默认值为0.23
    {
        m_dAET = 0;
        double dRLong, dRShort;
        dRLong = dRShort = 0;

        double dPrist, dPM;
        dPrist = dPM = 0;
        if (HydroSimulate.g_ModelRunPara.PETMethod == PET_METHOD.PET_PENMAN_MONTEITH)
        {
            PETInPristley prist = new PETInPristley(m_dTav, m_dHeight, m_nYear, m_nDn, 33);
            dRLong = prist.NetLongWaveRadiationRHmd(HydroSimulate.g_ClimatePara.slrdata.m_gridLayer.Values[row, col],
               HydroSimulate.g_ClimatePara.hmddata.m_gridLayer.Values[row, col]);
            dRShort = prist.NetShortWaveRadiation(dalbedo, HydroSimulate.g_ClimatePara.slrdata.m_gridLayer.Values[row, col]);
            dPrist = prist.PETByPristley(1.26, 0);
            m_dAET = 2 * dPrist - m_dPET;
            m_dPETPristley = dPrist;
        }
        else if (HydroSimulate.g_ModelRunPara.PETMethod == PET_METHOD.PET_PRISTLEY_TAYLOR)
        {
            PETInPM pm = new PETInPM(m_dTav, m_dHeight, m_nYear, m_nDn, 33);
            dRLong = pm.NetLongWaveRadiationRHmd(HydroSimulate.g_ClimatePara.slrdata.m_gridLayer.Values[row, col],
               HydroSimulate.g_ClimatePara.hmddata.m_gridLayer.Values[row, col]);
            dRShort = pm.NetShortWaveRadiation(dalbedo, HydroSimulate.g_ClimatePara.slrdata.m_gridLayer.Values[row, col]);
            dPM = pm.PETInPMByRHmd(HydroSimulate.g_ClimatePara.WindData.m_gridLayer.Values[row, col],
               HydroSimulate.g_ClimatePara.hmddata.m_gridLayer.Values[row, col], 0);
            m_dAET = 2 * m_dPET - dPM;
            m_dPETPristley = m_dPET;
        }
        else
        {
            PETInPristley prist = new PETInPristley(m_dTav, m_dHeight, m_nYear, m_nDn, 33);
            dRLong = prist.NetLongWaveRadiationRHmd(HydroSimulate.g_ClimatePara.slrdata.m_gridLayer.Values[row, col],
               HydroSimulate.g_ClimatePara.hmddata.m_gridLayer.Values[row, col]);
            dRShort = prist.NetShortWaveRadiation(dalbedo, HydroSimulate.g_ClimatePara.slrdata.m_gridLayer.Values[row, col]);
            dPrist = prist.PETByPristley(1.26, 0);

            PETInPM pm = new PETInPM(m_dTav, m_dHeight, m_nYear, m_nDn, 33);
            dRLong = pm.NetLongWaveRadiationRHmd(HydroSimulate.g_ClimatePara.slrdata.m_gridLayer.Values[row, col],
                HydroSimulate.g_ClimatePara.hmddata.m_gridLayer.Values[row, col]);
            dRShort = pm.NetShortWaveRadiation(dalbedo, HydroSimulate.g_ClimatePara.slrdata.m_gridLayer.Values[row, col]);
            dPM = pm.PETInPMByRHmd(HydroSimulate.g_ClimatePara.WindData.m_gridLayer.Values[row, col],
               HydroSimulate.g_ClimatePara.hmddata.m_gridLayer.Values[row, col], 0);
            m_dAET = 2 * dPrist - dPM;
            m_dPETPristley = dPrist;
        }

        return m_dAET;
    }
    /// <summary>
    /// 计算逐日叶面积指数
    /// </summary>
    /// <param name="method"></param>
    /// <returns></returns>
    double DailyLAI(DailyLAICalMethod method,GridLayerPara g_GridLayerPara)
    { 
        double dLai = 0.5;
        int VegID = (int)g_GridLayerPara.g_VegLayer.Values[row, col];
        if (method == DailyLAICalMethod.DAILY_LAI_CAL_SINE)
        {
            double dmax,dmin,doffset;
			dmax = g_GridLayerPara.VegTypes[VegID].LAIMX;
            dmin = g_GridLayerPara.VegTypes[VegID].LAIMN;
			doffset = g_GridLayerPara.VegTypes[VegID].doffset;
			dLai = ( dmax + dmin ) / 2.0 -( dmax - dmin ) *Math.Sin( 2 *Math.PI * ( m_nDn - doffset ) / 365) / 2;
        }
        else if (method == DailyLAICalMethod.DAILY_LAI_CAL_LINEAR)
        {
            throw new Exception("这种方法还没有实现");
        }
        else if (method == DailyLAICalMethod.DAILY_LAI_BY_MONTH)
        {
            dLai = g_GridLayerPara.VegTypes[VegID].LAI[m_nMon - 1];
        }
        return dLai;
    }
    /// <summary>
    /// 计算栅格逐日植被反照率
    /// </summary>
    /// <returns></returns>
    double DailyAlbedo(GridLayerPara g_GridLayerPara)
    {
        double dAlb = 0.23;
        int VegID = (int)g_GridLayerPara.g_VegLayer.Values[row, col];
        dAlb = g_GridLayerPara.VegTypes[VegID].Albedo[m_nMon - 1];
        return dAlb;
    }
    /// <summary>
    /// 计算栅格逐日盖度
    /// </summary>
    /// <returns></returns>
    double DailyCoverDeg(GridLayerPara g_GridLayerPara)
    {
        int VegID = (int)g_GridLayerPara.g_VegLayer.Values[row, col];
        return g_GridLayerPara.VegTypes[VegID].CoverDeg[m_nMon - 1];
    }
    /// <summary>
    /// 计算超过土壤田间持水量的土壤水下渗
    /// </summary>
    /// <param name="sw">土壤实际含水量</param>
    /// <param name="fc">田间持水量</param>
    /// <param name="dt"></param>
    /// <param name="dtPerco"></param>
    /// <returns></returns>
    double SWPercolation(double sw, double fc, double dt, double dtPerco)
    {
        double dret = 0;
        double swexcess = 0;
        if (sw > fc)
            swexcess = sw - fc;
        else
            swexcess = 0;
        dret = swexcess * (1 -Math.Exp(-1 * (dt / dtPerco)));
        return dret;
    }

    /// <summary>
    /// 计算无降雨时正常的时段土壤退水量
    /// </summary>
    /// <param name="sw"></param>
    /// <param name="dt"></param>
    /// <param name="dInfilRate"></param>
    /// <param name="dInitInfRate"></param>
    /// <returns></returns>
    double SWRecharge(double sw, double dt, double dInfilRate, double dInitInfRate)
    {
        double dret = 0;
        if (dInfilRate < dInitInfRate)
        {
            dret = dInfilRate * dt;
        }
        return dret;
    }														
	

	//栅格上各水文过程分量计算结果
    /// <summary>
    /// 地表径流
    /// </summary>
	public double	m_dSurfQ;	
	/// <summary>
	/// 壤中流
	/// </summary>		
	public double	m_dLateralQ;	
	/// <summary>
	/// 地下径流，基流
	/// </summary>	
	public double	m_dBaseQ;	
	/// <summary>
	/// 栅格总径流
	/// </summary>		
	public double	m_dTotalQ;				
    /// <summary>
    /// 冠层截流
    /// </summary>
	public double	m_dCrownInterc;			
	public double	m_dPET;	
	/// <summary>
	/// 实际蒸散发
	/// </summary>			
	public double	m_dAET;					
	public double	m_dAI;	
	/// <summary>
	/// 净雨量
	/// </summary>			
	public double	m_dNetRain;		
	/// <summary>
	/// 植被截留量赤字
	/// </summary>	
	public double	m_dCIDeficit;			
	
	public double	m_dInfil;				
	public double	m_dTbzero;				
	

	CanopyStorage	m_CanopyStore;	
	public double	m_dPETPristley;	
	
    /// <summary>
    /// 栅格降雨量值
    /// </summary>
	public double	m_dPcp;	
	/// <summary>
	/// 平均气温
	/// </summary>		
	public double	m_dTav;	
	/// <summary>
	/// 当前栅格列号
	/// </summary>		
	public int		col;	
	/// <summary>
	/// 当前栅格行号
	/// </summary>		
	public int		row;	
	/// <summary>
	/// 栅格的高程
	/// </summary>		
	public double	m_dHeight;			
	public int		m_nYear;				
	public int		m_nDn;	
	/// <summary>
	/// 土壤剖面稳定下渗率
	/// </summary>			
	public double	m_dFc;	
	/// <summary>
	/// 时段土壤下渗率
	/// </summary>			
	public double	m_dFp;	
	/// <summary>
	/// 当前处理时段的小时数
	/// </summary>
			
	public double	m_dHr;					
	/// <summary>
	/// 降雨强度
	/// </summary>
	public double	m_dRIntensity;			

	public int		m_nMon;				
	public double	m_dGridLAI;			
	public double	m_dGridCovDeg;
    /// <summary>
    /// 用于存放当前栅格中土壤实际含水量的变量，对应于SoilInfo类的SP_Sw变量
    /// </summary>
    public double Cell_SP_Sw;
};
}
